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Abstract: The eigenvalue distribution of Hoppe's two matrix model is investigated in 
detail as a function of the model's coupling. For small couplings it is a perturbed Wigner 



semicircle, while for large couplings it is a parabolic distribution which crosses over to 
a Wigner semicircle for eigenvalues within approximatley an inverse coupling from the 
boundary of the distribution. The model is approximately commuting at large couplings 
J> and we find the joint eigenvalue distribution of the two matrices. We also study a related 

three matrix model finding the corresponding three dimensional eigenvalue distribution 
there also. The techniques developed here are more widely applicable to other multi-matrix 
*• . models. 
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1 Introduction 

Multi-matrix models play an important role in several branches of modern physics especially 
in matrix string theory [1], the IKKT model [2] (and its lower dimensional variants [3]) 
and the BFFS and BMN models [4, 5]. They also describe the low energy dynamics of 
D-branes [6] and provide simple models of emergent geometry [7, 8] and emergent gravity 
[9, 10]. 

There are very few exactly solvable interacting multi-matrix models aside from Hoppe's 
two matrix model [11] which we analyse in this paper. This model plays a role similar to 
that of the two dimensional Ising model in critical phenomena, being exactly solvable for 
many quantities while having a rich behaviour that is not easily amenable to exact analysis. 



It also gives the characteristic behaviour one can expect in a wider class of multi-matrix 
models. The model was introduced as a guide to the physics of quantised membranes by 
Hoppe [11] and it has subsequently arisen in the low energy dynamics of D-branes [6] and 
in discussions of emergent geometry [12]. 

In [13] we studied the large coupling behaviour of Hoppe's model. We established that 
the matrices are approximately commuting at large coupling and the eigenvalue distribution 
of one of the matrices obeys a parabolic distribution. In an appendix we outlined the 
leading large coupling corrections to this parabolic distribution. We also established that 
the rotationally invariant 3-matrix model which reduces to Hoppe's model when one of the 
matrices is integrated out, has, at large coupling, a uniform joint eigenvalue distribution 



i 



within a ball of radius -R — ( ff ) where g is the coupling of the model. 

In this paper we pursue a more thorough investigation of Hoppe's model and make 
some further observations on its 3-matrix relative. 

The principal results of this paper are: 

• The eigenvalue distribution as a function of coupling, g, with perturbative expressions 
for large and small couplings. 

• The eigenvalue density of one matrix can be "lifted" to give rotationally invariant 
two and three dimensional distributions, which for large g become the eigenvalue 
distributions for the two matrix model and its three matrix relative, respectively. 

• The unique rotationally invariant "lift" to two dimensions is, for large coupling, a 
hemispherical distribution with a finite eigenvalue density at the boundary. 

• The unique rotationally invariant "lift" to three dimensions is the uniform distri- 
bution [13], however for any finite coupling, the "lifted" distribution grows in the 
shell 1/g ~ r < R and diverges at the boundary. When reduced to the one dimen- 
sional distribution this corresponds to a crossover to the Wigner distribution as the 
boundary is approached. 

The structure of the paper is as follows: 

In section 2 we introduce Hoppe's two matrix model, expand around a background of 
diagonal matrices, gauge fix so that one linear combination of the matrices is diagonal. We 
then integrate out the all remaining modes to obtain an effective action for these longitu- 
dinal modes (eigenvalues of the diagonalised matrix) . We then average over the arbitrary 
unit vector selecting the diagonalised matrix to get an integral equation for a rotation- 
ally invariant two dimensional distribution. We call this distribution the two dimensional 
"lift" of the eigenvalue distribution. The remainder of the section deals with analysing this 
"lifted" integral equation for weak and strong coupling. We establish that at weak cou- 
pling the solution is a uniform distribution while at strong coupling it gives a hemispherical 
distribution. The hemispherical lifted distribution in turn implies a parabolic eigenvalue 
distribution for the eigenvalue distribution of a single matrix. 

Section 3 considers the one dimensional eigenvalue distribution and begins by devel- 
oping perturbation theory around weak coupling. At zero coupling the eigenvalue distri- 



bution is the Wigner semicircle and we show that up to order (Rg) e (where R is the extent 
of the distribution and g the coupling) the distribution is a Wigner semicircle modified 
by polynomials in r\ = x/R (see equation (3.43)). Section 3.2 then develops perturba- 
tion theory for large coupling where the leading form of the distribution at large g is a 
parabola. An analytic form for the leading correction to the parabolic distribution is ob- 
tained. It is then shown that this reproduces the exact asymptotic growth of the observable 
v = g 2 < jf(X 2 ) > ,for large g, as obtained from and exact expression found in [6]. 

In section 3.3 we develop a numerical technique based on the Multhopp-Kalandiya 
method [14] to find the eigenvalue distribution for arbitrary couplings and verify our ana- 
lytic approximations for weak and strong couplings giving the regime of validity of these. 

Section 4 demonstrates that given a d — 1 dimensional distribution one can determine 
the rotational distribution "lifted" to d dimensions. It establishes that the uniform dis- 
tribution is the "lift" of the Wigner semi-circle and that lifting the parabolic distribution 
with its leading correction, equation (3.43), leads to a truncated hemispherical distribution 
(see figure 5). 

Section 5 discusses a 3 matrix variant of Hoppe's two matrix model which at large 
coupling was shown, [13], to have a eigenvalues uniformly distribution within a solid ball of 
radius (f^) ■ The integral equation for the rotationally invariant 3 dimensional "lifted" 
distribution is established. It is shown that quite generally the rotationally invariant d — 2 
dimensional distribution lifts to a rotationally invariant d dimensional distribution given 
by Pd(x) = — d 2^ x ■ The Wigner semicircle and perturbations of it at small couplings 
lift to distributions that are divergent at the boundary, while the parabolic distribution 
lifts to the uniform distribution within a ball. The leading corrections to the parabola lifts 
to a distribution that is again divergent at the boundary and we establish by numerical 
integration of the integral equation that the characteristic behaviour for all finite couplings 
is the square root divergence characteristic of the "lifted" Wigner semicircle. 

The paper has two technical appendices, the first, Appendix A, deals with the solution 
to the 2-dimensional "lifted" integral equation at large coupling while the second, Appendix 
B, deals with the large g asymptotics of the exact radial extent of the eigenvalue distribution 
R(g) and the large g asymptotics of the observable v. 

2 The two-matrix model. 

The principal model that we focus on in these notes is the two dimensional mass regulated 
model first considered by Hoppe [11]: 

Z= [ VX VYe- N ^ xl2+x22 -^ xl > x2 ^ . (2.1) 



Our main interest is the properties of this model at strong coupling, when it is in a nearly 
commuting 1 phase [12]. The main strategy to solve the model is to reduce it to a one 
dimensional model, from which many properties can be extracted exactly. 



^or large g we have < jf(i[X, Y]) 2 >~ ^ _ 5^(§f ) 2/3 > and so x and Y commute for g ->• oo. 



However we will first use the approach of ref. [13] and study the two matrix model 
directly by obtaining a two-dimensional distribution, which in the commuting phase co- 
incides with the joint eigenvalue distribution of the matrices. To this end we split the 
matrices as: 



-&-n — X. On + O,, 



Xi 



%i dij + CLij ; X i — [Xi , Xj ) ; (Xij — \ a ij i a ij ) 



Consider a constant unit vector n = (n , n 2 ) and define: 

x" = n(n.x); x = (1 — nn).x; a" = n(n.a); a = (1 — nn).a . 
Now we can use the SU(N) symmetry of the matrix model to fix the gauge: 

n.Sij = . 



(2.2) 



(2.3) 



(2.4) 



After integrating out the perpendicular elements of the matrices x and a the resulting 
effective action for (n.x) is [13]: 



1 N 1 N 

Ses[x] = — ^(n.fi) 2 -^^log 



8=1 



j,i=i 



?.X\2 



(n.(xi - Xj)) 



1 + g 2 {n.(xi - Xj)) 2 



(2.5) 



Next, we consider, for large N, the continuous limit of equation (2.5) and define a rota- 
tionally invariant two dimensional distribution p(x) so that (2.5) becomes 

1 f f j2 j2 / r->\ /-/m (n.(x - x')) 2 

- I I d z xd z x n(x)n(x') loer — v 



S eS [p(x)} = / d 2 xp(x)(n.x) 2 - 
+ fi [ I d xp(x) — 1 



1 + g 2 (n.(x — x')) 



+ (2.6) 



Note that in general the matrices X^ do not commute and p(x) is a rotationally invariant 
"lifted" version of the one dimensional distribution of n.x. However n.x is the eigenvalue 
of the matrix n.X and when the matrices commute diagonalising n.X would diagonalise 
both X 1 and X 2 . Therefore for g 2 — > oo when the model is in a commuting phase p(x) 
approaches the joint eigenvalue distribution of X^. At weak coupling the model is non- 
commuting and the "lifted" two dimensional distribution is a rotationally invariant lift of 
the eigenvalue distribution of one of the matrices, and is not itself an eigenvalue distribution. 
Varying with respect p in equation (2.6) we obtain: 



^ ^2 



u + (ft.x) 



d x p(x ) log 



(n.(x — x')) 



•/\\2 



1 + g 2 (n.(x — x')) 2 



(2.7) 



Note that equation (2.7) should be valid for any choice of n thus in order to obtain a 
rotationally invariant integral equation we average over n with weight one. 2 The result is: 



^ + 



d x p(x') ln(- 



x — X 



--)■ 



(2i 



1 + Vl + fF 



r>\2 



2 For ft = (cos (j>, sin <p) we integrate both sides of the equation by ^- J 

o 



2.1 2D distribution at weak coupling 



To obtain an integral equation suitable for perturbative calculation at small g we apply 
V 2 , on both sides of equation (2.6). We obtain: 



d x p(x 



4irp(x) - 2g 2 



V, 



^/i + g 2 {x 



P\2 



\X — X 
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-b{x — x) 



1 + a/1 + g 2 \x — x 



r'|2 



d x 



p(x>) 



(l+g 2 (x-x') 2 f/ 2 ' 



(2.9) 



which is an integral equation of the second kind. To avoid complications with the boundary 
of the integral (since the radius of the distribution runs with g) it is convenient to introduce 
a new variable if = x/R. Equation (2.9) can then be written as: 



p(v) 



1 (%) 2 

2vr 2vr 



j2 / 
d 7] 



P(ff) 



W\<i 
The kernel in equation (2.10) is bounded by: 

K 2 (Rg,\ri-r?\) = ~ 



(i + (R g y(ff-ff 



\2^3/2 ' 



f] 



(l + {Rgf^-ffff 2 



< 1 



(2.10) 



(2.11) 



therefore for sufficiently small (Rg) 2 we can solve equation (2.10) iteratively: 



P(jj) 



1 

2^ 



1 + (Rg) 2 i d 2 r]'K 2 + (Rg) 4 i d 2 ^K 2 I d 2 tf'K 2 + . . . 



(2.12) 



Using equation (2.12) one can solve for p(rj) perturbatively to arbitrary order in (Rg). Here 
we present the solution to sixth order: 

Note that R in equation (2.13) is g dependent and expansion in g would look different. We 
can determine the g dependence of R using the definition of i] and the normalization of 
p(x) = p(rj). Indeed: 

d 2 xp(x) = R 2 d 2 f]p(r]) = 1 (2.14) 

We can now substitute equation (2.13) in equation (2.14), expand R in terms of g 2 and 
determine the coefficients in the expansion by solving equation (2.14) order by order. This 
procedure can be performed to arbitrary order. Here we present the result to sixth order 
in g: 



r 1 2 15 4 165 
R=V2 -T2 9+ W=2 9 "8V1 



9 6 + 0(g s ) . 



(2.15) 



2.2 2D distribution at strong coupling 

Next we focus on the large g limit of the model. It is convenient to modify equation (2.8) 
to: 

// + ^ = 2 /dVp(*')ln( f~ A ). , (2.16) 

J \ + \j\ + g 2 \x-x'\ 2 

where the factor of g on the right-hand side of the equation is compensated by a redefinition 
of the constant p — > p! . It is a straightforward exercise to obtain the asymptotic form of 
equation (2.16) in the g — > oo limit: 

M' + f = --fd 2 x> J ^\- + 0{l/g 2 ) . (2.17) 

2 g J \x — x'\ 

To leading order the integral equation that we obtain is: 

li + - = ~ 2 [Mt0&- = - 2 Lv 7&f> , P{X>) _ , (2.18) 

2 g J \x-x'\ g J J J x i + x '2 _ 2xx' cos 6 

o o v 

where we have used p(x) = p(\x\) = p{x). After performing the integral over (j) we arrive 
at the integral equation: 

v ' 

where K(z) is the complete elliptic integral of the first kind. The integral equation (2.19) 
can be solved [15] (see also Appendix A) for p(x): 



9 \{R 2 -p')-x* 








o 


^ VR 2 - x 2 





p{x) = ^ ^_ - = JL ^R2 _ X 2 (2 . 2 o) 

where we have fixed the constant p! = —R 2 by demanding that the distribution be finite 
at the boundary (x = R). Equation (2.20) is the hemisphere distribution of ref. [12]. By 
normalizing the distribution p(x) to one we can fix the radius of the distribution it!: 

3^ 1/3 



R =[ -) . (2.21) 

Having obtained the hemisphere distribution for g — > oo directly in two dimensions we are 
interested in the behaviour of the model for finite values of the coupling constant g, when 
the model is nearly commuting. Note that strictly speaking the correction to the joint 
eigenvalue distribution of the model at finite coupling is not well defined since the model 
is not in a commuting phase. However, there is a complex observable $ = X 1 + iX 2 , 
which has complex eigenvalues which are well defined at any coupling. Furthermore in 
the commuting phase the real and imaginary components of the eigenvalues of $ coincide 
with the components of the joint eigenvalues x. Numerically one can simulate the model 
at finite g and obtain the distribution of <3? keeping in mind that in the commuting phase 



this is the joint eigenvalue distribution. In order to compare to numerical simulations we 
need to understand the behaviour of the distribution for large but finite g. 

It turns out that it is technically easier to determine the correction to the hemisphere 
distribution by integrating out one of the matrices and study the corresponding one di- 
mensional distribution. In the next section we analyze the reduced model and the one 
dimensional distribution at general coupling g. As we show by solving an integral equation 
of Abel's type we can lift the one dimensional distribution to a rotationally invariant two 
dimensional one. 

3 The one matrix model 

In this section we focus on the one-dimensional distribution of the matrix n.X defined in 
the previous section. Without loss of generality we can choose n = (1,0) and x = (x,y). 
The integral equation for the distribution pi(x) is given by: 



// + x 2 = / dx' pi(x') log 



g 2 (x - x') 2 
1 + g 2 (x — x') 2 



(3.1) 



where we have substituted the definition of pi(x): 

VR 2 -x 2 



Pi(x)= / p(y/x 2 + y 2 )dy . (3.2) 



-V-R 2 -x 2 

in equation (2.7) and have redefined the constant p — >• p' to add the factor of g 2 in the 
argument of the logarithmic function in equation (3.1). It is convenient to differentiate 
equation (3.1) with respect to x. The resulting integral equation can be written as: 

R R 

dx' ^^ + f dx' Pl (x') K(g, x'-x) , (3.3) 



-R -R 

where the kernel K(g, u) is given by: 

K(g,u) = --^r^. (3.4) 

1 + g z u z 

3.1 ID distribution at weak coupling 

At g = we have K(0, u) = and the integral equation (3.3) has a simple Cauchy kernel: 

-x=ldx' P 4^. (3.5) 

J x — X 

-R 

Note that at vanishing coupling the model is Gaussian and hence the one dimensional 
distribution of X 1 , or equivalently n.X, should be a Wigner semi-circle. The unique 
bounded solution of equation (3.5) is indeed the Wigner semicircle (3.6). 

px (x) = -\/R 2 -x 2 , (3.6) 

7T 



with radius R = v2 ; which agrees with the Rg —> limit of equation (2.15). The pertur- 
bative solution for small g can be obtained, in terms of the coordinate rj = x/R and the 
coupling Rg, by integrating out one of the components of ff in equation (2.13) : 



p(m) = / p{\] fii + ril) dm , 

where p(rj) is related to pi{x) via: 

p(v) = -j=pi(Rri) ■ 



(3.7) 



(3.8) 



The final expression for small Rg up to sixth order is given by: 
Pi(v) = \/l-?7 2 



1 (Rg) 2 (4t? 2 + 3) (Rg) 4 (8t? 4 + 20r? 2 + 9) (Rgf 
ir 2ir 8ir 16ir 



+ O ((Rg) 8 ) 



(3.9) 

As one can see, the small Rg corrections deform the Wigner semicircle but it still has the 
characteristic \J 1 — i] 2 behaviour. 

3.2 ID distribution at strong coupling 

For large g one can use: 



to obtain: 



l - + K(g,u) = --5'(u) + 0(l/g 2 ) 
u g 



x=- Pl (x) + 0(l/g 2 ) , 



which to leading order in g is solved by the parabolic distribution [12]: 



(3.10) 
(3.11) 

(3.12) 



with radius given by equation (2.21). 

In order to obtain the corrections to the parabolic distribution at finite g we will 
derive an integral equation of the second kind which can (at least in principal) be solved 
iteratively. Let us begin by noting that the kernel of the integral equation (3.3) can be 
written as: 

(3.13) 



— h K(g, x - x) = —Ki(g,x'-x) , 

x — x ax 



where K\(g, x' — x) is the symmetric kernel: 



Ki(g,x' -x) = -log 



(3.14) 



g 2 (x — x') 2 
1 + g 2 (x — x') 2 

After integration by parts, and noting that p(±R) = 0, equation (3.3) can be written as: 

R 
x= dx' Ki(g, x — x) p'\(x) . (3.15) 

-R 



Again, to deal with the g dependence of the limits of the integral, it is convenient to rewrite 
the integral equation (3.15) in terms of the variables rj = x/R, Rg and the distribution 
p(rj) = ^ so that (3.15) becomes: 

i 
ri = J difK x (Rg,rf-ri)pHjt). (3.16) 

-l 

At large Rg the kernel K\ has the expansion: 

Kx = ~8{jf - V ) + {l/{Rgf) . (3.17) 

Kg 

Next we define the kernel: 

AK(Rg, r/ - r?) = -^-K x {Rg, n ' - rj) - S(rf - rj) (3.18) 

7T 



with Ap defined by: 



where pry, given by: 



Ap'{rj) = p'{rj) -RgP(o)(v) = p(v) + — V , (3-19) 

Po(ri) = ^(l-V 2 ), (3-20) 

is the parabolic distribution (3.12) valid in the strict g — >• oo (Rg —¥ oo) limit. The integral 
equation (3.16) can be written as: 

l l 

Ap'(ri) = -Rg J drf AK(Rg, rj - V ) p[ 0) (r/) - j dr,' AK(Rg, V ' - r,) Ap'(r,') . (3.21) 

-l -l 

Equation (3.21) is an integral equation of the second kind for the correction Ap. Further- 
more, from the definition of AK and the expansion of K\ at large Rg, equation (3.17), 
it follows that AK dies out at large Rg and hence the integral equation (3.21) can be 
developed recursively in a convergent series: 

Ap'( V ) = -Rg J AK p{ 0) + Rg j AK j AK p{ 0) - Rg J AK J AK J AK p[ 0) + ... . 

(3.22) 
If we define pLJRg,r)) = —Rg f AK p, Q s and: 

i 
p{ n+1) (Rg,rj) = - fd V AK(Rg,ri'- V )p{ n) (Rg, V ') n = l,2,..., (3.23) 

-l 

we arrive at the following expression for Ap(ri): 

A*^E/*/%^ ) . (3--) 



71=1 i 



At large Rg we have 1/Rg < J AK < 1. Therefore: 

1 



(Rg) 



n-l 



P{n)( R 9,V) 



(Rg) 



n—l 



(Rg) 



71-2 



(3.25) 



and naively one would expect that at large Rg the contribution to Ap in equation (3.24) 
from terms with n > 1 would die out. However one can show that p',Joo,rj) is not 
integrable near the boundary (rj = ±1) and by regulating it with a cutoff of the order 
~ l/(Rg) one can estimate that for large Rg: § pt,ARg,rj) ~ [Rg) n for n > 1 and 
J /5'qs (Rg, 7/) ~ log(-R<j) for n = 1. Therefore for n > 1 we have that p', n JRg,r])/(Rg) n ^ 1 ~ 
[(5(1 — rj) — 5(1 + 77)] for sufficiently large -Rg and all terms with n > 1 give a constant 
contribution n ~ 1 to Ap in equation (3.24) as long as 77 £ (— 1, 1)- At the boundary the 
contribution from all terms vanishes and we have Ap(±l) = 0. Therefore to leading order 
we have the following expression for Ap: 



Ap(ri) 



J d V 'p'(Rg,r,') + « + 0((\og(Rg)/Rg) if - 1 < 7/ < 1 
-1 

if 7/ = ±1 



(3.26) 



Using the definitions from equations (3.20) and (3.23) we can obtain the following expres- 
sion for fjLy 



p' {1) (Rg,v) = ^r [tan" 1 [Rg(l + V )]+ tan" 1 [Rg(l - r,)} - vr] + ^ lo S 



- ^ (1 - r/ 2 ) log 



4vr 2 



Note that: 



P(i)(°°^) 



;i + T ? ) 2 (i + (i t > g ) 2 (i-T ? ) 2 ) 
;i-7 / ) 2 (i + ( j r 5 ) 2 (i + t ? ) 2 ) 

1 — 7/ 



1 7/ 1 

"vr 2 l-7/ 2 + 2^ 2 g 



1 + 7/ 



l + (flg) 2 (l->7) 2 - 

l + ( J R 5 ) 2 (l + 7/) 2 

(3.27) 



(3.28) 



which is indeed not integrable near 77 = ±1. Let us introduce a cutoff e. The regulated 
expression for pm is: 



pffifa) 



1 

2^ 



7/ log 



1 — 7/ 

1 + 7/ 



+ 2^ l0g 



+ O (e log(e)) . 



(3.29) 



The cutoff e can be expressed in terms of Rg. Indeed if we integrate directly the integrable 
expression (3.27) we obtain: 



p {1) (Rg,rj) = ^[T 1 (Rg,r ] )+T 2 (Rg, V )+T 3 (Rg,r ] )} , 
where T\, T2 and T3 are given by: 

tan-^Rgtl - 7?)) + tan^i^l + 7/)) 



(3.30) 



T 1 (Rg,7 1 ) = (l- V 2 ) 1 



7T 



tan- 1 (2.Rg) - tan~ 1 (^ ff (l + 7/)) - tan~ 1 (^(l - y)) ,_ „. 
+ 37r(ifo)2 ' (3 - 31) 
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T 2 (Rg,v) 



2vr Rg 



r\ log 



l + (ifr) 2 (l-T/) 2 

1 + (Rg) 2 (l + r/) 2 



+ log(l + 4( J R 5 ) 2 ) 



(3.32) 



13(^,77) = |f log 



(^) 4 (l-r? 



2\2 



(l + ( J R 5 ) 2 (l-r ? ) 2 )(l + (^) 2 (l + r ? ) 2 ) 



Rg 

2tt 



V 1 



log 



2M 



(l-r ? ) 2 (l + (fig) 2 (l+7?) 

(l+7 ? ) 2 (l + ( J R 5 ) 2 (l- ?? ) 2 ) 

2flj/,__r 4(i? 5 ) s 



3tt 



log 



1 + 4(i^ 



For the large Rg expansion of pn\ we obtain: 



P(i)(Rg,v) 



2vr 2 



77 log 



1 — 77 

1 + ??. 



+ ^ log(2 e 3 / 2 %) + 0(log(%)/(%)) 



(3.33) 



(3.34) 



Comparing equations (3.29) and (3.34) we conclude that e = e~ 3 / 2 / (Rg) ~ l/(Rg), which 
agrees with the analysis performed below equation (3.25). Taking into account the constant 
k from equation (3.26) we arrive at the following expression for p in the interval 77 € (—1,1): 

1 



p(v) 



Rg 

2vr 



y 1 ~ V ) + SZ2 V log 



2^ 2 



V 



I + 77 



+ ^ log(2 e 3/2 Rg) + « + 0(\og(Rg)/(Rg)) . 



(3.35) 

Note that at large i?g the correction, A/5, vanishes at 77 = 1 — <5i, where to leading order 
<5i = — 2TrRn n an< ^ W (z) is the Lambert's product log function -- the solution to 
z = W(z)e ( z > . Therefore, at large Rg, to a very good approximation Ap(ry) is given by 
equation (3.35) in the interval 77 G (—1 + <5i, 1 — Si) and can be taken as zero outside this 
interval. 

The constant k in equation (3.35) can be determined numerically by computing pr n ) 
up to sufficiently large n. However one can also fix k indirectly by comparing to some of 
the exact relations for the observables of this model. Indeed, following the approach of 
ref. [6] one can obtain an exact relation between the radius of the distribution and the 
coupling constant g. This relation is rather complex and is given parametrically in terms 
of elliptic integrals. For a more detailed derivation we refer the reader to appendix B. Here 
we provide only the first few terms in the large g expansion of R. 



R 



(3tt 



1/3 



g 



_ 1/3 _ 2l0gg + l0g(967T 4 ) g _ 1 + q/ -5/3) 

6lT 



(3.36) 



The radius of the distribution R can also be calculated using that: 



1 



drjp(rj) 



Rg 

2ir 



di](l 



v 2 ) + 



drj Ap(r/) . 



-1 



-1 



-1 



(3.37) 



Substituting the expression for A/5 from equation (3.35) into equation (3.37) and solving 
for R order by order in g one arrives at: 



R 



(*L\ 1/a n -v* 



V 2 J 



g 



41og5 + 21og(127r) + 12vr 2 K + 3 _ x 



12vr 



g-' + Oig- 5 ' 3 ) 



(3.38) 
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Comparing equations (3.36) and (3.38) we obtain: 

log(4vr 2 ) - 1 



4tt 2 



and we find the correction Ap is given by: 



&p(n) = 7^2 ^ lo s 



1 — 77 

1 + ?? 



log(47rflg) + l / togCggA 

27T 2 +ty V ^ y 



(3.39) 



(3.40) 



and our expression for the p(rj) in the large Rg regime is 



p(v) 







1 - v 2 ) + 2^ V lo S 



1=2 



+ 



log(47rRg)+l 
2tt 2 



if |7?| > ] - S 



+ OfS^gfiij if |r?| < 1 



(3.41) 
where 5 = W(l/e)/(27rRg) ~ l/i?g and W(.z) is the Lambert's product log function. 

Our expression for the leading correction to the distribution can be tested by calcu- 
lating the observable v = g 1 < TrX 2 >, which was obtained in closed form in ref. [6] (see 
also appendix B). Indeed from the definition of v and p{r\) it follows that: 



R 



v = g 2 f dxx 2 pi(x) = R 2 (Rg) 2 f d^piv) 



(12vr) 2 / 3 



,4/3 



20 



-R 



(12vr) 2 / 3 



9 2/3 + 0(g c 



(3.42) 

which is in perfect agreement with the result of ref. [6]. In deriving (3.42) we have used 
equations (3.36) and (3.41). 

Finally let is obtain the leading order behavior of pi(x) at large g. Using p\{x) = 
Rp(x/R) and equations (3.36) and (3.41) we obtain: 



Pi{x) 



9_ 

2vr 



(if- 



X 

+ 2^ l0g 



3vr\ 

igj 



1/3 



3tt\ 

2g) 



1/3 



+ X 



+ i(S) 1/3+o (T |(3 - 43) 



In the next subsection we develop a numerical routine to obtain the one dimensional 
distribution for arbitrary values of the coupling (Rg). 

3.3 Interpolating solution for general coupling. 

In this subsection we construct an interpolating solution to the integral equation: 



1 

/ dr] ' „j _„ + / dr]/ pM) K (Rg, v' - v) 



rj' — rj 



(3.44) 



where K is defined in equation (3.4). Equation (3.44) is a singular integral equation of the 
first kind with a Cauchy kernel. An interpolating solution of this equation can be found 
using the Multhopp-Kalandiya method (see ref. [14] Chapter 14.5). When the source is odd 
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(—rj in our case) the approximate solution bounded at the boundary of the distribution 
(77 = ±1) is given by: 



p a (cOS0) 



vr(2n- 



tt£* 



costy j 



^(ra 1 SUI 



nlir 
2n+l 



sin((2n + 1)0) 



1=1 



cos(2#) — cos 



2ln 
2n+l 



(3.45) 



where cos# = 77, 2n is the number of nodes into which the interval [—1, 1] is divided and Q\ 
are related to the roots of the Chebyshev polynomial of the second kind: 



rji = cos 61, 61 



In 



1 = l,...,2ra . 



(3.46) 



2n + l ' 
The value of p at the nodes rji is determined by the system of linear algebraic equations: 



^2(ck,i ~ C2n+l-k,l)p( cos ^l) 

1=1 



cosOk, k = 1, . . . n , 



(3.47) 



where 



Ckl 



sin #2 
2n+l 



2efcz 



cos 9i — cos 6k 



+ K(Rg, cos 9k, cos ( 



and eik = (k - I) mod 2 . (3.48) 



By choosing sufficiently large n one can generate numerical solution of almost arbitrary 

precision. Using: 

1 

-1 



R 



dr]p(r]) 



and pi(x) = Rp{x/R) 



(3.49) 



one can generate an approximate numerical solution for p\{x). In figure 1 we present the 
plot of the distribution p\{x) for range of coupling constants < g < 133.3. The black 
dashed curve in the figure represents the Wigner semi-circle (3.6) one can see the perfect 
fit with the data for g = 0. The red dashed curve fits the curve for the highest value of 
g (g = 133.3) and represents the analytic expression for the distribution at large g given 
in equation (3.43). One can see the excellent agreement of the numerical results with the 
analytic analysis from the previous subsection. 

Let us now verify the approximate solutions for p(rj) obtained in the previous section. 
In figure 2 we have presented plots of p{rf) for small values of the coupling constant < 
Rg < 0.6. The dashed red curves in figure 2 represent the approximate solution for small 
Rg from equation (3.9) while the blue ones are the numerical solution. One can see that 
the approximation is excellent for Rg < 0.6 and is reasonably good for Rg = 0.6. 

Next we focus on the large Rg regime. In figure 3 we present plots of p(rj) for 3 < 
Rg < 18. The dashed red curves represent p(rj) given in equation (3.41). One can see how 
the approximation improves as we increase Rg and at Rg = 18 it is already excellent. 

To further verify the correctness of our numerical approach we evaluate some of the 
observables of the model which can be obtained in closed form [6] (see also appendix B). 
In figure 4 we present plots of the radius of the distribution and the observable v (defined 
in equation (B.2)) as a function of the coupling constant g. 
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Figure 1. The eigenvalue density p\(x) for small coupling < g < 133.3. The black dashed curve 
is the Wigner semi-circle (g = 0) while the red dashed curve represents the approximate solution 
(3.43) for g = 133.3 i.e. Rg ~ 42.745. 
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Figure 2. The figure shows the eigenvalue distribution pi(x) for < g < 0.6. The dashed red curves 
represent the approximate solution (3.9) while the blue curve represents the numerical solution for 
the corresponding value of g. The bottom curve is the Wigner semi-circle g = 0. 

The blue curves in figure 4 represent the exact results for the radius R and the observ- 
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3 < Rg < 18 




Figure 3. The figure shows p{rj) in the large Rg regime < Rg < 18, with the red dashed lined 
the theoretical expression (3.41) and the blue curves the numerical solution. 




Figure 4. R(g) and v(g) The blue curves are the exact results (B.5) and (B.2). The dashed 
curves are the large g approximate expressions (B.8) and (B.4), while the doted red curves are the 
numerical solution to (3.44). 

able v given by equations (B.5) and (B.2). The dashed curves are the large g approximate 
expressions (B.8) and (B.4) one can see that for g > 1 there is excellent agreement with 
the exact result. Finally the doted red curves represent the numerical result obtained by 
solving numerically the integral equation (3.44). One can observe the perfect agreement of 
the numerical and the exact results. 

4 2D distribution at general coupling 

In this section we present our numerical results for the lifted 2D distribution at general 
coupling. The easiest way to achieve this is to "lift" the distributions of the previous 
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section. We will describe this in some generality. 

4.1 Lifting the distribution and Abel's integral equation 

Let us comment on the relation between rotationally invariant distributions in different 
dimension, which is a generalization of the relation (3.2). Consider a d— 1-dimensional dis- 
tribution pd-i obtained by reducing a d-dimensional rotationally invariant distribution p^: 

VR 2 ~r 2 

Pd-i{r)= I p d (Vr 2 + z 2 )dz . (4.1) 

-VR 2 -r 2 

Using the change of variables: r = \/R 2 — ( , z = \/£ — r/ the integral equation can be 
written as: 

«. l( ^0-/«^*|. (4-2) 



Equation (4.2) is Abel's integral equation with a solution given by [? ]: 

*(„) _ I - / O^jE^M -±'"f Sptt * . (4.3) 

rK ' irdriJ s/rf^t v =R 2 -w 2 irooduj ^/ r 2 - uj 2 

o R 

It is an easy exercise to obtain the two dimensional distribution corresponding to p\ via 
equation (4.3). At zero coupling the Wigner semicircle (3.6) corresponds to the a uniform 
distribution p = ^-, while at strong coupling the parabolic distribution (3.12) reproduces 
the hemisphere distribution (2.20). We will use equation (4.3) to "lift" the approximate 
solution to p\ for general couplings. 

4.2 Numerical results 

Using equation (4.3) we can "lift" the numerical solution from the previous section for the 
one dimensional distribution to obtain a two-dimensional rotationally invariant distribution 
which in the commuting phase of the model coincides with the joint eigenvalue distribution. 
Furthermore for large coupling, g, we can lift the approximate expression for p\ from 
equation (3.43). Note that since \x\ < R ~ g -1 ' 3 the second term in equation (3.43) is of 
order log(g)/g 1 ' 3 and dies out at large g. The lift of the first term is part of a semicircle 
of radius (37r/2) 1 ' 3 5 _1 ' 3 . The lifted distribution at large g is then given by: 



p(x) 



\ 2/3 \ 1 / 2 / s 

^1 -x 2 ) +0{ l m\ , for0<z<ii 



, for x > R 



where R is given by equation (3.36). Note that at the boundary (x = R) the distribution 
is non-zero. Using equation (3.36) one can estimate the magnitude of the distribution at 
the boundary: 

(R) = >A)S(9faV ) i/» + (45) 

VW 2 1 /63l/3 7r 7/3 y ^ ■ • • ' ^ °> 
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which is growing with g. However the magnitude of the maximum of the distribution, at 
x = 0, is p(0) = (37r 5 /2) 1 ' 5 g 2 ' 3 which grows faster with g. Therefore in the limit g — > oo 
the shape of the distribution approaches a hemisphere. 

In figure 5 we present our numerical results for the distribution p(x) for different 
coupling constants < g < 12. The continuous curves represent our numerical results for 




1.2 



1.4 



Figure 5. The figure shows the lifted two dimensional distribution as the coupling is increased. 
The dashed black curve (uniform distribution) is the lift of the Wigner distribution at g = 0. The 
red dashed curves are the approximation (4.4) while the solid curves are the numerical results for 
the lifted distribution. 

the lifted distribution and the colour changes from bright blue to violet as g increases. The 
black dashed curve represents a uniform distribution of magnitude l/(2ir), one can see the 
perfect fit with the numerical results for g = (the light blue curve). The red dashed 
curves represent the approximate expression (4.4) and one can see how the approximation 
improves as g increases and at g = 12 (Rg ~ 8) it is already excellent. 

5 Three matrix model realization 

Let us now consider the model (also originally introduced parenthetically by Hoppe [11] 
page 73): 

Z = f VX VYVZe- N ^ x2+Y2+z2 - ia ^ Y ^ . (5.1) 

It is easy to verify that if one integrates out the Z matrix and defines g 2 = (ia) 2 /A 
one recovers the two matrix model (2.1). This suggests that the model (5.1) should be 
as solvable as the two matrix model. Note also that there is a global SO (3) symmetry 
rotating the X, Y and Z matrices. In ref. [13] the model (5.1) was analyzed in the spirit 
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of section 1 . Namely the Hermitian matrices X, Y and Z were split into diagonal and off 
diagonal modes: 

ij — '^i"ij'Q'iji *ij — ^i^ij'^ij] ^ij — Xi"ij~T~Q'ij> %i — [%% > %% > X )] Q>ij — \®ij ' ®ij ' ®ij ) 

(5.2) 
and introduced an axial gauge n.a = 0, where n is a three dimensional unit vector. After 
integrating out the perpendicular degrees of freedom one arrives at the following effective 
action: 



1 N 
S eS [(n.x)} = j-f^2(n-X; 



i=l 



1 

2N 2 



N 



E 1 ^ 



ij=l 



g 2 {n.(xj - Xj)) 2 
1 + g 2 (n.(xi - Xj)) 2 



+ ^_l) log92 . (5 .3, 



Then considering a coarse grained approximation and varying the corresponding distribu- 
tion function p 3 one arrives at the equation: 

g 2 (n.(x - x")) 2 



^2 



p + (ft.x) 



d A x' '/93(f) log 



1 + g 2 (n.(x — x')) 



(5.4) 



Equation (5.4) can be averaged over a unit two-sphere (integrating both sides of the equa- 
tion by jr J ^2) to obtain: 



!^2 



d 3 x' p 3 (x*) 



-2arctan(g|x — x'| 



+ log 



g 2 (x — x') 2 
1 + g 2 (x — x') 2 



g\x — x I 
Next we apply the Laplacian A x to both side of equation (5.5). The result is: 

P3(X') 



1 



3 / 



d 6 x 



(5.5) 



(5.6) 



\x — x'\ 2 {\ + g 2 \x — x'\ 2 ) 

where x = \x\ and we have used that the lifted distribution p 3 is rotationally invariant 
(p 3 (x) = p 3 (x)). To obtain p 3 we need to solve the integral equation (5.6). Integrating 
over the angular coordinates in (5.6) and multiplying both sides of the equation by x results 

in: 

R 

(x + x') 2 (l + g 2 (x - x') 2 ) 



I dx ir x p 3 {x') 



log 



(x-x') 2 {l + g 2 {x + x') 2 ) 



(5.7) 



If we extend the integral in equation (5.7) over an even interval the integral equation can 
be written as: 



dx' (-2-kx 1 'p 3 (\x'\)) K x (g,x - x) , 



(5.8) 



-ft 



where K\ is the kernel (3.14) from section 3. Comparing equations (5.8) and (3.15) one 
arrives at the following relation between p\ and p 3 : 



P3(x) 



Pi 0*0 ? 

2lTX 



x > . 



(5.9) 



In fact equation (5.9) can be proven in more generality. Let us consider a rotationally 
invariant distribution in d dimensions. The rotationally invariant distribution in d — 2 
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dimensions obtained by integrating out two of the spacial dimensions can be obtained by 
integrating over a disk: 

VR 2 -x 2 R 

p d _ 2 (x) = 2ir j Pd (Vx 2 + r2)rdr = 2TT f ^ d( , (5.10) 

x 

where in the last expression we defined C, = Vx 2 + r 2 . Now after differentiating the first 
and the last expressions in (5.10) by x, the integral equation for p d reduces to an algebraic 
one which can easily be solved to obtain the analogue of equation (5.9): 

Pd{x) = -%^ , x > . (5.11) 

These considerations confirm that the procedure of averaging over n in equation (5.4) is 
equivalent to lifting the one dimensional distribution (via equation (5.11)). 

It is a straightforward exercise to "lift" the results of section 3 to the three dimensional 
case using equation (5.9): 

5.1 3D distribution at weak coupling 

At vanishing coupling the Wigner semicircle (3.6) is lifted to: 

P3W = i vrb • (5J2) 

which is divergent but integrable at the boundary. In analogy with the one dimensional 
case where the distribution behaves as V R 2 — x 2 near the boundary for any finite cou- 
pling, equation (5.12) suggests that the lifted three dimensional distribution will diverge 
as 1/V R 2 — x 2 for any finite coupling g. 

To obtain the perturbative expression for p% at small g it is convenient to change 
variables r\ = x/R and /^(r?) = Rps(Rr]). The lift of equation (3.9) is then given by: 



fe(v) 



VT^ 



1 (Rg) 2 (Urf - 5) (Rg) 4 (40?? 4 + 28r? 2 - 31) {Rgf 

(5.13) 



p + ! g- ''"" iJ'"" + ' !L>"""" +0 « fl "> 8 > 



5.2 3D distribution at strong coupling 

In the limit g — > oo the one dimensional distribution is parabolic (3.12) with radius (2.21). 
The lifted three dimensional distribution is uniform (obtained in ref. [13]): 

f%( X ) = ^2 OT foil) = ^2 • ( 5 - 14 ) 

Note that since the model is commuting in the limit 5—7-00, the lifted distribution in 
equation (5.14) is also the three dimensional eigenvalue distribution of X, Y, Z. 
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It is straightforward to lift the correction to the distribution Ap (3.41) at large but 
finite g: 



fc(j)) 




2vr 3 (l-r ? 2 ) 



4^ ^g 



1=3 



+ o(mm) if |i,i < i 

if M > 1 



<5 

5 



(5.15) 



with (5 = W(l/e)/(2nRg) as in (3.41). Finally we lift the approximate expression for p\ 
(3.43) to obtain the expression for p%: 



P3(x) 



g t J_ 
2vr 2 + 2tt 3 



1/3 



2 9 J X 



47T 3 : 



log 



1/3 



1/3 



+ X 



+ o 



logfl 
5 



(5.16) 



where x E (0, i?) and R is given in equation (2.21). 

5.3 3D distribution at general coupling 

In this subsection we use equation (5.9) to lift the interpolating solution p from section 3. 
In figure 5 we present a plot of the numerical solution for the expression y 1 — rj 2 p^ (77) 
for Rg e [5,55]. The color of the curves changes from blue to violet as Rg increases. The 
red dashed curves represent the approximate expression (5.15) for y 1 — rj 2 Ps(tJ). One can 
see that the approximation improves as one moves far from the boundary 77 = 1. One can 
also observe how the approximation improves as Rg grows. 



p(W 



1-77 2 



3.0 




Figure 6. The figure shows the three dimensional rotationally invariant lifted distribution ^3(77) 
rescaled by \f\—r\* for 5 < Rg < 55. The continuous curves are the numerical solutions to the 
integral equation and the dashed red curves are the approximate solution (5.15). 
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6 Discussion 

We have performed a rather detailed study of Hoppe's 2-matrix model (2.1) and in particu- 
lar developed both perturbative and interpolating solutions for the eigenvalue distribution 
of either matrix. We "lifted" this one matrix eigenvalue distribution, i.e. the one dimen- 
sional distribution, to rotationally invariant distributions in both two and three dimensions. 
For large couplings these lifted distributions capture the joint eigenvalue distributions of 
the two and three matrix models. 

We found that the two dimensional distribution does not go to zero at the boundary but 
rather has a finite value, p(R), given in (4.5), that grows as p(R) ~ <?s ln(g) for large g; see 
Figure 6. This implies that the distribution lifted to a rotationally invariant 3-dimensional 
distribution must diverge at the boundary for any finite g. 

From Figure 6, and the fact that this is the lift of the two dimensional distribution 
shown in Figure 5, we can deduce that near the boundary the asymptotic behaviour of 
Psix) is given by 

»(x) ~ £- 2 + . MR ! . (6.1) 

The divergence 3 of p%(x) as x approaches the boundary is an essential feature as without it 
the limiting two dimensional distribution could not attain a non-zero value at its boundary. 
This in turn means that the one dimensional distribution, asymptotically close to the 
boundary, is given by 

Pl {x) ~ 2p(R)y / R 2 -x 2 (6.2) 

and the distribution crosses over to a Wigner semicircle as the boundary is approached. 

The implication of this is that the noncommutative modes, which of necessity are 
present in the model, are concentrated near the boundary of the distribution, i.e. they are 
associated with the largest eigenvalues. 

The techniques used in this paper are applicable to a wide variety of models. The large 
coupling analysis of section 3, taking advantage of the <5-convergent sequence nature of the 
kernel, is novel and is easily applicable to a wider class of models. Also the interpolation 
approach of section 3.3, is novel and can be adopted to more general situations. 

The models studied here have the important feature that for large couplings the dom- 
inant configurations are commuting matrices. The fluctuations around these commuting 
modes are always present and so the full matrices never truly commute. This is in part 
reflected in the divergence, as the boundary is approached, of the "lifted" three dimensional 
distribution. This divergence is subleading for large coupling but an essential feature of 
such models none the less. 

3 In contrast to the approximation (5.16), which is not integrable at the boundary, the form (6.1) is 
integrable with the divergence capturing the nature of the distribution in the region 1 — 5 < J^ < 1. 
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A Solving for the hemisphere distribution 

In this appendix we solve the integral equation (2.19). Substituting y(x) = xp{x) we arrive 
at: 

«*>~§("4)-/<$H!^)- (A ' 1) 

x / 

The solution to equation (A.l) is given by [14, 15]: 

*> « < 7iM* m= i(°m±. ( A.2) 



vr 2 dx J Vt 2 - x 2 ' dt J V^-s 



— o 2 

o 



Substituting f{x) from equation (A.l) into equation (A. 2) we obtain: 

y(x) = x — ft — 2 , (A. 3) 

which implies equation (2.20). 

B Exact Results 

In this appendix we provide with slight extension some of the exact results for the model 
(2.1) obtained in ref. [6]. One of the exact results of the authors of ref. [6] was a closed 
form expression for the observable: 

R 

v = g 2 [ dx Pl (x) . (B.l) 

-R 



It is given by: 

1 K 2 10tf 2 (i? + m-2) + 2#(6-6m + m 2 ) + (l-m)(m-2) 
U{m '~12~5^ M 2 + 2(m - 2)i9 + 1 - m ' ( ' 

where K = K(m) and •& = E(m)/K{m) [E and K are the standard elliptic integrals). The 
elliptic modulus m can be determined in terms of the coupling constant g via: 

g 2 (m) = — (_3# 2 + 2(2-m)tf- (1 - m)) . (B.3) 

37T 4 
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Equations (B.2) and (B.3) specify (in parametric form) the g dependence of the observable 
v. For large g one can obtain the expansion: 

_ (127T) 2 / 3 4/3 3 2/3 / 1 1 \ 2/3 

Another exact result obtained in ref. [6] relevant to our discussion is the radius of the 
distribution R for which the authors derived the following integral presentation: 

R=- f dt X3 ~ f , (B.5) 

X4 

where xi, X2, X3 and X4 are functions of m given by: 

Xl = A-(2-m-20); Z2 = A- (1 - 2i?); (B.6) 

5 Z 7T Z g( Z 7T Z 

^ 2 K 2 

x 3 = -s— ^ (3$ + m — 2); x 4 = -~— ^ (1 - m - 2i?); . 
5 z 7r z 5 7T Z 

The integral (B.5) can be solved in closed form: 



„. . K(m) , 1 1 — •&{m) , ,„ . 

R(m) = \ \ Z sin" 1 \ \-L\ m ) (B.7 

it gym) V m 

where Z(<j) \ m) is the standard Jacobi Zeta function. To obtain a large g expansion of R 
we expand equation (B.7) near m = 1. Using equation (B.3) we can obtain the expansion 
used in equation (3.36): 

^ (| y/ % _ 1/3 _ 2 ,o g9 + M 9fa% . 1 + _^ 9 . 5/3+0(9 _ 7/3) (B8) 
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